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We extend a lattice Boltzmann algorithm of liquid crystal hydrodynamics to include an applied 
electric field. The approach solves the equations of motion written in terms of a tensor order param- 
eter. Back-flow effects and the hydrodynamics of topological defects are included. We investigate 
some of the dynamics relevant to liquid crystal devices; in particular defect-mediated motion of 
domain walls relevant to the nucleation of states useful in pi-cells. An anisotropy in the domain wall 
velocity is seen because defects of different topology couple differently to the flow field. 
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I. INTRODUCTION 



The coupling of the optic and electric response of liquid crystals has lead to their wide application in display devices 
in recent years. Liquid crystalline materials are often made up of long, thin, rod-like molecules This anisotropy 
is what leads to their useful optic properties. The molecular geometry and interactions can lead to a wide range of 
equilibrium phases. In this paper we are concerned with the nematic phase where the molecules tend to align along a 
preferred direction, referred to as the director, giving long-range orientational order. The propensity to order, as well 
as the direction along which the system orders are conveniently described by a tensor order parameter Q . 

An important feature of the nematic phase is the possibility of topological defects. These are director configurations 
that cannot relax to the nematic ground state by a continuous rotation of the molecules. Examples of such defects 
that we will consider are shown in Fig. ^(b). In addition, although widely recognized as an important factor in device 
performance, the flow of the liquid crystal in response to changes in the orientation of the molecules is both difficult 
to measure experimentally and to incorporate into a simulation of the device. 

Liquid crystal hydrodynamics are typically described by the Ericksen-Leslie-Parodi equations of motion [Q . However 
these are restricted to an order parameter field of constant magnitude and do not include the hydrodynamics of 
topological defects. Here we consider the Beris-Edwards formulation of liquid crystal hydrodynamics This allows 
for variations in the magnitude of the nematic order parameter and correctly models both defect dynamics and the 
coupling between the velocity field and the motion of the order parameter. 

The Beris-Edwards equations are complex and previous numerical work to explore different fiow regimes is very 
limited. In Ref. these equations were used but the problem was simplified by imposing a velocity field and ignoring 
back-flow (back- flow refers to the effect of the order parameter dynamics on the velocity field) . In Q an Euler solution 
of a somewhat simpler version of the Beris-Edwards model was used to study the effect of hydrodynamics on phase 
ordering in liquid crystals. However it has recently been shown that the Beris-Edwards equations can be solved 
successfully using a lattice Boltzmann algorithm Q . 

In this paper we will extend the lattice Boltzmann approach described in Q to include the electric field terms 
necessary for a description of electro-optic device physics. The model is also generalized to include the three elastic 
constant formula for the deformation free energy of the liquid crystal. The equations of motion are presented in 
Section II and an outline of the numerical algorithm in Section III. 

We then apply the algorithm to some examples of interest in device physics ||^. The nematic liquid crystal is 
confined between two plates a few /im apart. The director configuration on the plates is fixed. When an electric field 
is switched on, the molecules align parallel or perpendicular to the field depending on the sign of the anisotropy of 
the dielectric constant. After switching off the field, long-range elastic interactions ensure that the molecules reorient 
themselves in the direction preferred by the surfaces. The device can be used as a display because different liquid 
crystal orientations have different optical properties. 

Switching between different director orientations can take place homogeneously throughout the bulk of the device 
or it may proceed through domain growth involving the motion of domain walls, together with associated defects 
The operational state of a device such as a pi-cell may be topologically distinct from its state at zero voltage. Before 
the device can be used, the operational state must be nucleated and must grow to fill the display. It is important 
to examine how this depends on different parameters of the liquid crystal material. In particular liquid crystals are 
complex fluids and the speed of domain motion will depend on their hydrodynamic properties. 



1 



In section IV we examine the global stability of two topologically distinct states as a function of the applied voltage 
and the tilt angle at the plates. In Section V we discuss the relevance of defect motion to domain growth between 
the two states and show that hydrodynamic coupling significantly changes the defect velocities. 



II. THE HYDRODYNAMIC EQUATIONS OF MOTION 

We summarize the formulation of liquid crystal hydrodynamics described by Beris and Edwards Q and extended to 
include electric fields. The continuum equations of motion are written in terms of a tensor order parameter Q which 
is related to the direction of individual molecules m by Qa0 = {mamp — ^5oip) where the angular brackets denote a 
coarse-grained average. (Greek indices will be used to represent Cartesian components of vectors and tensors and the 
usual summation over repeated indices will be assumed.) Q is a traceless symmetric tensor. Its largest eigenvalue, 
|g, < g < 1, describes the magnitude of the order. 

We first write down a Landau free energy which describes the equilibrium properties of the liquid crystal. This 
appears in the equation of motion of the order parameter, which includes a Cahn-Hilliard-like term through which 
the system evolves towards thermodynamic equilibrium. 

Free energy: The equilibrium properties of a nematic liquid crystal can be described by the Landau-de Gennes 
free energy 

= / dV {fbulk + f elastic + f field] + I dS {fsurf} ■ (1) 

Jv JS 
fbulk describes the bulk free energy 

hulk = -^(1 - ■|)Qq/3 J^QatiQfSfQ-fa H ^iQaf3)'^- (2) 

For 7 = 2.7 there is a first-order transition from the isotropic to the nematic phase, f elastic is the analogue of the 
Frank elastic free energy density 

/elastic = -^{daQpj)'^ + {daQ af){d0Q j3j) + -^Q a/sidaQ je){d(3Qje) ■ (3) 

The dielectric constant is usually measured along and perpendicular to the nematic axis, n, and is usually assumed 
to give a relation between the electric displacement D and field E of the form 

D = eiE+(e|| -ei)(n-E)rl (4) 

For a uniaxial nematic, this is entirely equivalent to assuming that 

2 

ea/3 = gCaQa/S + (5) 

where 

3 

ea = ^(e|| -Ci), (6) 

e™ = ge± + g^ll (7) 

The electric contribution to the thermodynamic potential f field is 

If 1 1 

f field — / D • C?E = — — e,„i? — Y^f-aEaEpQap. (8) 

At the surfaces of the device we assume a pinning potential 

fsurf = -^O^siQaP — Qap)^- (9) 

This corresponds to a director at the surface preferring to lie along the direction of the eigenvector of Q*^ corresponding 
to the largest eigenvalue 2/3g. 
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Equation of motion of the nematic order parameter: The equation of motion for the nematic order parameter 

is i 

(^t+u- V)Q-S(W,Q) = rH (10) 

where F is a collective rotational diffusion constant. The first term on the left-hand side of equation ( |l0| ) is the 
material derivative describing the usual time dependence of a quantity adverted by a fluid with velocity u. This is 
generalized by a second term 

S(W, Q) (CA + f2)(Q + 1/3) + (Q + I/3)(eA - ^) 

-2^(Q + I/3)Tr(QW) (11) 

where A — (W + W-^)/2 and O = (W — 'W^)/2 are the symmetric part and the anti-symmetric part respectively of 
the velocity gradient tensor Wap = dpUa- S(W, Q) appears in the equation of motion because the order parameter 
distribution can be both rotated and stretched by flow gradients. This is the consequence of the rod-like geometry of 
the liquid crystal molecules. ^ is a constant which will depend on the molecular details of a given liquid crystal. 

The term on the right-hand side of equation ( [lo| ) describes the relaxation of the order parameter towards the 
minimum of the free energy. The molecular field H which provides the driving motion is related to the derivative of 
the free energy by 

= Hbulk + Helastic + Hflgld + Hgurf , (12) 



where 



Hbulk = -^o(l - ^)Q + Ao7 (Q^ - (I/3)rrQ2) - Ao7Q^rQ^ (13) 



Helastic,a.f3 — 

1 , r„ „ „ . „ w„ „ . 1 



1277^"*^ 3 

Hsurf = -as(Q-Q"), (16) 



Hf,eid,»p ^ ^{Ea^Ep - ^E^^), (15) 



and Sai3 is the Kronecker delta. For (Q) dz — is assumed, and the symmetry and tracelessness of Q is also exploited 
for simplification. 

When computing the functional derivative in ( p^ ) surface terms arise from the integration by parts. For reference, 
these are 

^^s«r/ ^ / dsSQafi {Liajid^Qafi) + '^L2[aaid^Qyi3) + <Ti3{dyQ^a:)] + Lsa^Q^^id^Qa/s)} , (17) 



dV 



where a is the surface (unit) normal. 

Continuity and Navier-Stokes equations: The fluid momentum obeys the continuity 

dtp + dapua = 0, (18) 

and the Navier-Stokes equation 

pTf 

pdtUa + pUpdpUa = dpTap + OpCTap + —dp{{5ap - 3dpPo)djUj + d^Up + dpUa), (19) 

where p is the fluid density and tj is related to the viscosity. The form of this equation is not dissimilar to that 
for a simple fluid. However the details of the stress tensor reflect the additional complications of liquid crystal 
hydrodynamics. There is a symmetric contribution 
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1 sj~ 

+ '2^{QaP + -:^Sai3)QjeH~^^ - dpQ^^ h aM,af3 (20) 

and an antisymmetric contribution 

Tq/S = Qa^H^p — Ha^Q^p. (21) 

The background pressure Pq is taken simply to be 

Po = pT, (22) 

and is found to be a constant in our simulations to a very good approximation. The Maxwell stress tensor (TM,af3 |lC| ] 
describes the stress due to the external electric field 

0'M,a8 = T;—{DaEp + EaDi3 ~ DjEjSap). (23) 



III. A LATTICE BOLTZMANN ALGORITHM FOR LIQUID CRYSTAL HYDRODYNAMICS 

We now describe a lattice Boltzmann algorithm which solves the hydrodynamic equations of motion of a liquid 
crystal (10), (|8|), and (|l9|). Lattice Boltzmann algorithms are defined in terms of a set of continuous variables, 



usefully termed partial distribution functions, which move on a lattice in discrete space and time 

The simplest lattice Boltzmann algorithm, which solves the Navier-Stokes equations of a simple fluid, is defined in 

terms of a single set of partial distribution functions which sum on each site to give the density. For liquid crystal 

hydrodynamics this must be supplemented by a second set, which are tensor variables, and which are related to the 

tensor order parameter Q ||] . 

We define two distribution functions, the scalars fi{x) and the symmetric traceless tensors Gi{x) on each lattice 

site X. Each fi, Gi is associated with a lattice vector e^. We choose a nine- velocity model on a square lattice with 

velocity vectors et = (±1,0), (0,±1), (±1,±1), (0,0). Physical variables are defined as moments of the distribution 

function 

P^J^J^^' PUc.^^^^e^a. Q = ^G,. (24) 

i % i 

The distribution functions evolve in a time step At according to 

j,(x + e,Ai, t + Ai) - /,(f , <) = ^ [C/.(x, + C/,(f + e,At, t + At, {/*})] , (25) 

(f + e; At , i + At) - Gj (f , t) = 

^ \CG^(x, t, {G,}) + CG^{x + e,At, t + At, {G*})] . (26) 

This represents free streaming with velocity and a collision step which allows the distribution to relax towards 
equilibrium. /* and G* are first order approximations to fi(x + CiAt, t + At) and Gi(x + eiAt, t + At) respectively. 
Discretizing in this way, which is similar to a predictor-corrector scheme, means that lattice viscosity terms are 
eliminated. In a lattice Boltzmann algorithm for a simple fluid these terms can be added to the physical viscosity 
term. Here this is not the case: the lattice viscosity gives an additional, spurious term in the equations of motion. 
An additional advantage of the predictor-corrector approach is that the stability of the lattice Boltzmann algorithm 
is improved. The collision operators are taken to have the form of a single relaxation time Boltzmann equation [Q, 
together with a forcing term 

C/,(f,t,{/,}) - -l(/,(f,t) -/;«(x,t, {/,})) +p,(f,t, {/,}), (27) 

Cg.(x, t, {G,;}) = --(G,(f , t) - G^'?(f , t, {G,})) -f M,(x, t, {G,}). (28) 
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The form of the equations of motion and thermodynamic equihbrium follow from the choice of the moments of the 
equilibrium distributions f^"^ and G^^ and the driving terms pi and M^. f^"^ is constrained by 

i i i 

where the zeroth and first moments are chosen to impose conservation of mass and momentum. The second moment 
of f'^'' controls the symmetric part of the stress tensor, whereas the moments of pt 

'^Pi=0, ^PiCia = df^Tafl, ^PiCiaCif^ = (30) 
i i i 

impose the antisymmetric part of the stress tensor. For the equilibrium of the order parameter distribution wc choose 

G^' = Q, ^ Gl'^eia = Qm„, ^ Gl'^e^aeip = QwaW/?- (31) 

i i i 

This ensures that the order parameter is convected with the flow. Finally the evolution of the order parameter is 
most conveniently modeled by choosing 

^M, = rH(Q) + S(W,Q) = H, ;^M,e„ = (^M,)«a, (32) 

i i i 

which ensures that the fluid minimizes its free energy at equilibrium. 

Conditions (|29|)~(|3^) can be satisfied as is usual in lattice Boltzmann schemes by writing the equilibrium distribution 
functions and forcing terms as polynomial expansions in the velocity |^ 

ft'' = As + BsUaCia + C gU^ + D sUaUpeiaCifj + EsapB^aCiiS, 
G^' = Js + KsUaCia + LgU^ + N sUaU/seiaet/S , 
Pi — TgdfjTapeia, 

Mi = Rs + SsUQCiQ., (33) 

where s = Ci ^ g {0, 1, 2} identifies separate coefficients for different absolute values of the velocities. A suitable choice 
is 



A2 = 


-{(yxx + 


(Tyy)/16, Al = 


2.42, Ao^p-l2A2, 


B2 = 




Bi = 4^2, 




C2 = 


-p/16, 


Ci = -p/8, 


Co = -3p/4, 


D2 = 




Di = p/2 




E2xx 


~ i'^vv ~ 


0'a;a;)/16, E2yy 


= ^E2xx, E2xy = E2yx 


E\xx 


= 'iE2xx, 






Jo = 


Q 






K2 = 


= Q/12, 


Ki = 4K2, 




L2 = 


-Q/16, 


Li = -Q/8, 


Lo = -3Q/4, 


N2 = 


= Q/8, 


Ni = Q/2 




1^2- 


1/12, 


Ti = 4r2, 




R2 - 


:H/9, 


Ri = Rq = R2 




S2 = 


H/12, 


Si -4S2, 





(34) 

where any coefficients not listed are zero. 

IV. H-V PHASE DIAGRAM 

We now use the lattice Boltzmann algorithm described in Sections 2 and 3 to investigate a simple liquid crystal 
device. The device consists of a nematic liquid crystal confined between two planes a distance apart. At the 
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surface Q*^ is chosen to give a surface tilt to the director of +6p at a: = and —9p at x = L^- At zero apphed voltage 
these conditions result in a global minimum free energy state with a splayed director configuration, or H state as 
shown in Figure 0(a). At high voltages, typically on order of 6V, the H state is no longer the global minimum, and 
a bend configuration (V state) is obtained like the one shown in Figure 0(b). At intermediate voltages, the V state is 
more relaxed, like the one shown in Figure ^(c) . 

The V state may remain for some time even at zero voltage as it is metastable. As the H and V states are 
topologically distinct, the transition from one to the other requires nucleation and the generation of defects. This 
process is illustrated in Figure |2[ 

The simulations were performed on a 48 x 48 grid. We used Li = 0.0440, L2 = 0.0445, L3 = 0.0606, ea = 41.4, 
em = 9.8, Aq = 1.0 and 7 = 3.0. By a choice of appropriate length, time and pressure scales these correspond to 
physical values values Li = 17 A pN, L2 — 17.6 pN, L3 = 24.2 pN, = Ly = 3 /im which is consistent with the 
nematic liquid crystal E7 (Merck Ltd. UK) at 25°C. Periodic boundary conditions were used in the y direction and 
bounce-back in the x direction. 

A phase diagram for the H and V states in the tilt angle/voltage plane is shown in Fig. ^. The system was started 
in a configuration like that shown Figure 0. For low fields the H domain has a lower free energy density, thus it grows 
(crosses). For high enough fields the H domain shrinks (circles). The boundary is defined as the voltage at which 
the two domains have the same free energy. This diagram compares well with those measured experimentally for a 
system using E7 jist . 



V. DOMAIN GROWTH 



Consider the dynamics of the domain growth in the system. Rather than fix parameters to a specific material, 
we wish to explore how the dynamics are affected by different parameters. In particular, we are interested in how 
hydrodynamics affects the speed of the domain growth. This is motivated by the observation in Ref. |l^ that the 
domain growth was anisotropic and the speculation that this may be due to hydrodynamics. To isolate effects due to 
different elastic constants, for these simulations we set L2 = = 0.0, the so-called one-elastic constant approximation. 

The initial configuration, depicted in Fig. ^(a), is a horizontal (i.e., along the y-direction) domain in an otherwise 
vertically aligned state. This models a time shortly after the electric field has been switched off when small but 
macroscopic domains have formed in the device. 

Once the simulation commences the director configuration relaxes rapidly to that shown in Fig. 2|(b). The horizontal 
and vertical domains are both tilted slightly by the elastic coupling to the surface spins. Then defects are formed at 
the domain boundaries, which merge to form single defects with strengths ±1/2 at the centre of each domain wall 
respectively. Once the two defects have formed the vertical domain begins to grow and the defects move in opposite 
directions. Our aim is to investigate the effect of back-fiow and surface director tilt on this motion. Most previous 
work on the dynamics of topological defects has ignored the back- flow [pT|-p^ . 

We first investigate the effect of the surface tilt dp on the defect speed. The horizontal domain grows because it 
decreases the free energy of the device. It can be shown |l^ that the difference in free energy between the horizontal 
and vertical domains is proportional to 45° — dp, where dp is the angle of the surface tilt to the y-axis. Thus as 9p 
increases, the free energy difference decreases, and the defects move more slowly. At 9p = 45° the two domains have 
the same free energy and the defects stop moving. For 9p > 45° the horizontal domain begins to shrink. 

A particular advantage of the simulation is that the back-flow can easily be switched off by setting CTq/j — —PoSap 
and Tai3 to zero in (^) and (^J). Thus the effect of the back- flow can be unambiguously identifled. Since no flow is 
imposed, this leads to a zero velocity field throughout the whole simulation. The dynamical equation for the system 
in this case can be obtained from (nw by setting u to zero. It corresponds to a simpler model which does not include 
hydrodynamics. 

Consider first the diamonds in Fig. ^. This corresponds to the case with back-flow switched off. For this case both 
defects move at the same speed (but in opposite directions). This can be explained by considering a local co-ordinate 
transformation x — > —x which is equivalent to changing the sign of the off-diagonal elements in Q. Assuming m = 0, 
the dynamical equation (|l^) is invariant under this transformation whereas the two defects in Fig. ||(b) transform 
into each other. 

The triangles and circles in Fig. show the velocity of the defects when back-flow is included in the model. Now 
the two defects move with different speeds since the xy elements of the stress tensors (^0|) and (^) are not invariant 
under the x — s- —x transformation. Thus the stress flelds are different for the two defects. Due to this stress field two 
different fiow fields are formed around the defects (see Fig. ||). The flow is stronger around the s = +^ defect and it 
is substantially accelerated. The speed of the s — defect is affected much less by the flow. 
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An experimental setup similar to Fig. g was considered in ||T^ where the growth of a circular twist domain in a 
horizontal environment under the influence of an external field was studied for different surface tilt values. It was 
found that the domain growth was anisotropic. Although this problem is three-dimensional, a vertical cross section 
through the growing domain gives a geometry similar to that considered here. It seems extremely plausible, that the 
essential physics is captured by our model with the anisotropy in growth resulting from the difference in the flow fields 
corresponding to different defect topologies at the defect boundary. Quantitative mappings between the experiment 
and the simulation results must be treated with caution because of the different dimensionality. However, in both the 
simulation and the experiment the difference between the wall speeds was about 40 — 50%. 

The simulations were performed on a 700 x 48 grid. We used Li — 1.0, L2 — — 0, Aq — 5.0, 7 = 3.0, F — 0.44, 
^ = 0.52 and Tf = 1.0. The larger Li corresponds to increasing the resolution of the lattice in order to perform 
accurate defect velocity measurements. 



VI. SUMMARY AND DISCUSSION 



In this paper we have described a lattice Boltzmann algorithm to simulate the dynamics of a non-Newtonian fluid, 
liquid crystals. 

In particular we have considered the the presence of an external electric field and the possibility of a general Frank 
free energy with three elastic constants. In the continuum limit we recover the Beris-Edwards formulation within 
which the liquid crystal equations of motion are written in terms of a tensor order parameter The equations are 
applicable to the isotropic, uniaxial nematic, and biaxial nematic phases. Working within the framework of a variable 
tensor order parameter it is possible to simulate variations in the magnitude of order and hence the dynamics of 
topological defects and non-equilibrium phase transitions between different flow regimes. 

The algorithm was used to explore states in nematic liquid crystal devices. We find that switching between 
topologically distinct states is, as expected, strongly inhibited by a free energy barrier between the initial and the 
final states. In this case the state can only change once a domain of the final state has been nucleated. Defects form 
at the moving domain walls and hence this geometry allows an investigation of defect hydrodynamics. We find that 
a s = +i defect is substantially speeded up by back- flow effects. This is relevant to device physics where the speed 
of switching is an important design variable. 

There are many directions for further research opened up by the rich physics inherent in liquid crystal hydrodynamics 
and the generality of the Beris-Edwards equations. For example, another possible switching pathway would be via 
the escape of the director into a third dimension. The addition of flexoelectric terms to the equations of motion will 
allow problems relevant to bistable liquid crystal displays to be addressed. Even cases where the switching is the 
consequence of back-flow [|l6| can be modeled. Numerical investigations are vital to address these problems because 
of the complexity of the equations of motion. 
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FIGURE CAPTIONS 

Fig. 1 Alignment of molecules for a tilt angle +9p on the top surface and —dp on the bottom surface. Director 
configuration when 9p < 45° and (a) the field is switched off and the system has had time to relax to its global 
minimum; (b) the field is switched on the a fairly high voltage ~ 6V; and (c) the field is at a low voltage ~ 2V or 
lower. The system may remain in the metastable state (c) for some time even at zero voltage. 

Fig 2. (a) Initial director configuration used to study domain growth. After the vertical electric field is switched off, 
a horizontal domain is created in order to nucleate the growth process, (b) Shortly thereafter two defects are formed 
at the boundary of the horizontal and vertical domains. The left (right) defect has a topological strength s = 
(s ~ +^)- Both the horizontal and the vertical domains are slightly distorted due to the defects and the surface tilt. 

Fig. 3. Phase diagram for the H and V states as a function of tilt angle and voltage. At low fields (crosses) the H 
domain has a lower free energy density, thus it grows. If the field strong enough (circles), the H domain shrinks. On 
the line the two domains have equal free energy density. 

Fig. 4. The velocity of the two defects as the function of surface tilt if back-flow is ignored (diamonds) and included. 
Notice that if back-flow is not included then the two defects move with the same speed. Hydrodynamics accelerates 
the s = +^ defect (triangles) substantially, while it affects less the s — defect (circles). 

Fig. 5. (a) Velocity field corresponding to the director configuration shown in Fig. g(b). Notice that there is a 
strong vortex structure around the s = +^ defect and at the defect core it points in the direction of defect movement. 
The flow at the core of the other defect is much weaker, (b) Flow field corresponding to the dashed box. For better 
visibility the length of the vectors is changed disproportionally. 
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